#pragma once

#include <cmath>
#include <limits>

#define SIGN(a) ((a >= 0) ? 1 : -1)
#define ITERATION 100

class Function
{
protected:
    double (*f)(double); // 函数
    Function *D;         // 函数链，存放函数

public:
    Function(double (*f)(double)) : f(f), D(nullptr) {}
    Function(double (*f)(double), double (*g)(double)) : f(f), D(new Function(g)) {}

    double operator()(double x) { return f(x); }                                      // 获取指定点函数值
    double derivation(double x) { return f(x + FLT_MIN) / FLT_MIN - f(x) / FLT_MIN; } // 获取指定点的导数值
    double next(double x)
    {
        if (D != nullptr)
        {
            return (*D)(x);
        }
        return 0;
    }

    double bisection(double u, double v);  // 二分法
    double Newton(double a, double b);     // 牛顿法
    double secant(double x_0, double x_1); // 割线法

    ~Function()
    {
        if (D != nullptr)
        {
            delete D;
            D = nullptr;
        }
    }
};

double Function::bisection(double u, double v)
{
    double u_f = f(u), v_f = f(v);
    if (SIGN(u_f) == SIGN(v_f))
    {
        return 0;
    }
    double mid, mid_f;
    // 限制最多迭代100次
    for (int i = 0; i < ITERATION; i++)
    {
        mid = (u + v) / 2;
        mid_f = f(mid);
        // 根据精度判断退出
        if (fabs(u - v) < FLT_MIN || fabs(u_f - v_f) < FLT_MIN)
        {
            break;
        }
        else if (SIGN(u_f) == SIGN(mid_f))
        {
            u = mid;
            u_f = f(u);
        }
        else
        {
            v = mid;
            v_f = f(v);
        }
    }
    return mid;
}

double Function::Newton(double a, double b)
{
    double dx, root = (b + a) / 2;
    for (int i = 0; i < ITERATION; i++)
    {
        if (next(root) == 0)
        {
            return root;
        }
        dx = f(root) / next(root);
        if (fabs(dx) < FLT_MIN)
        {
            break;
        }
        root -= dx;
    }
    return root;
}

double Function::secant(double x_0, double x_1)
{
    double dx;
    for (int i = 0; i < ITERATION; i++)
    {
        if (f(x_1) - f(x_0) == 0)
        {
            return x_1;
        }
        dx = f(x_1) * (x_1 - x_0) / (f(x_1) - f(x_0));
        if (fabs(dx) < FLT_MIN || fabs(x_1 - x_0) < FLT_MIN)
        {
            break;
        }
        double tmp = x_1;
        x_1 -= dx;
        x_0 = tmp;
    }
    return x_1;
}

#undef SIGN
#undef ITERATION